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ABSTRACT 

Astrophysical shocks or bursts from a photoionizing source can disturb the typical collisional plasma 
found in galactic interstellar media or the intergalactic medium. The spectrum emitted by this plasma 
contains diagnostics that have been used to determine the time since the disturbing event, although 
this determination becomes uncertain as the elements in the plasma return to ionization equilibrium. 
A general solution for the equilibrium timescale for each element arises from th e elegant eige nvec- 
tor method of solution to the problem of a non-equilibrium plasma described by [Masai ( 1983) and 
Hughes & Hclfand (1985). In general the ionization evolution of an element Z in a constant electron 
temperature plasma is given by a coupled set of Z+l first order differential equations. However, 
they can be recast as Z uncoupled first order differential equations using an eigenvector basis for the 
system. The solution is then Z separate exponential functions, with the time constants given by the 
eigenvalues of the rate matrix. The smallest of these eigenvalues gives the scale of slowest return to 
equilibrium independent of the initial conditions, while conversely the largest eigenvalue is the scale 
of the fastest change in the ion population. These results hold for an ionizing plasma, a recombining 
plasma, or even a plasma with random initial conditions, and will allow users of these diagnostics to 
determine directly if their best-fit result significantly limits the timescale since a disturbance or is so 
close to equilibrium as to include an arbitrarily-long time. 
Subject headings: atomic data, atomic processes, plasmas 
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1. INTRODUCTION 

Astrophysical plasma densities can have densities in 
the range of a few particles per cubic meter, resulting 
in interaction timescales of millions of years that create 
long delays between a thermodynamic event and its ap- 
parent consequences. For example, the ion populations in 
remnants created by the passage of supernova shocks are 
often well out of equ ilibrium with the electron temper- 
ature in the plasma (|Itohl Il978t iGronenschild fc Mewe I 
11983 Hamilton. Sarazin fc Chevalier I I1983D . The ge"rT 
eral problem of calculating the evolution of such a 
plasma as a function of time, including the emission 
spectrum from the plasma during the evolution, re- 
quires a large collection of atomic data and a moder- 
ately complex spe ctral code. A number of such codes 
exist, such as the iRavmond fc Smith I (fl 977) code, the 
iKaastra. Mewe fc Nieuwenhuiizen I (|1996t ) (SPEX) code, 
and thelSinitheLalJ ((20010 (APEC) code. 

Comparing observed X-ray spectra with these spec- 
tral models have been used to d etermine the 'age ' 
of supernova rem nants (SNR) (e.g. lHarrus et all Il997t 
lYatsu et al.1 l2005h or in the case of Cas A, to de- 
termine the when individual knots in the ejecta were 
shocked and thus follow the propaga tion of the origi- 
nal explosion ([Hwang fc Laming~ll2003| ). Although these 
examples are for ionizing plasmas, recent Suzaku ob- 
servations have also shown th at recombining plasmas 
can be found in S NR as well (lYamaguchi et al.l 120091 : 
IQzawa et al.l 120097 ). Fitting the spectra of these non- 
equilibrium plasma returns the temperature and overall 



density- weighted timescale of the plasma n e xt (because 
all ionization and recombination rates scale with electron 
density n e ), but these fits contain no direct indication of 
how close this is to an equilibrium value for each element. 
As a result, the spectroscopist may estimate a result that 
is in fact best treated as a lower limit. The differences 
in elemental timescales amount to 'clocks' that operate 
over various ranges as each element returns to equilib- 
rium, so with high-resolution spectra the sensitivity of 
the timescale measurement can be enhanced by focussing 
on specific elemental lines. 

While modifying the APEC code to use a faster and 
more general ionization state calculation routine (de- 
scribed in Ej2j) , we noted that one useful side effect of this 
method was a direct measure of the timescale for obtain- 
ing both the first changes to the ion population and the 
final equilibrium ionization state for each atom regardless 
of initial conditions and assuming only that the plasma 
is held at constant electron temperature. These results 
are shown graphically in and some consequences are 
given in 21 



2. METHOD 

IMasail (|l98l and lHuehes fc Helfand I (fl985l) described 
a fast method to calculate the ionization population of 
an atom using the eigenvectors and eigenvalues of the 
rate matrix equati on. Restating for convenie nce the for- 
malism set up in [Hughes & Hclfand (1985), the basic 
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equation is: 
dF- 

= n e (a i -i(T e )F i - 1 -[a i (T e )+R i - 1 (T e )}F i +R i (T e )F i+1 ) 

(!) 

where Fi is the ion fraction of the ith ion with i = 1 the 
neutral case and i = Z + 1 the fully-stripped ion with 
atomic number Z. All the rates are proportional to the 
electron density n e , which has been factored out. If n e is 
not constant, this approach is still applicable using n e t 
as the independent variable. The ionization rate ai(T e ), 
a function of the electron temperature T e , gives the rate 
from state i to state i + 1, while the recombination rate 
Ri(T e ) gives the rate from state i+ 1 to i. This equation 
can then be rewritten in matrix form as 

^ = -n e A(T e ).F (2) 

where F is a vector with dimension Z + 1 of ion popula- 
tions and A(T e ) is the matrix of ionization and recombi- 
nation rates. This matrix will have one zero eigenvalue, 
representing the equilibrium solution. This singularity 
may cause numerical problems that can be addressed us- 
ing the normalization requirement that ^ Fi — 1 to re- 
duce the matrix to Z x Z d imensions (see the Appendix of 
IHughes fc Helfand I (fl985h for details). Since the method 
assumes the electron temperature is held constant, it is 
convenient to define G = F — F eq , where F eq is the ion 
population in equilibrium at temperature T. Then we 
can rewrite Eq[2]as 

dC 

^ = -n e A(T e ) • G, (3) 

dropping the constant term A(T e ) • F eq as dF/dt = 
in equilibrium. This formulation has the advantage that 
lim G = by definition. 

t— foo 

Then we can define V(T e ) as the matrix of eigen- 
vectors of A(T e ) with eigenvalues A. Then setting 
G' = V(T e ) _1 G and A = AI, we can write 

with the solution G'(t) = cexp(— n e Xt) where c is deter- 
mined from the initial conditions. The values of G' do 
not, of course, correspond to individual ions, although 

these can be easily recovered using G = VG". Expand- 
ing this for ion i, we see that 

Gi(t) = ]T VjiG'j = Vim exp(-n e \ 3 t). (5) 
j j 

The requirement that lim G = implies that A 7 - > for 

i->oo 

all j. More importantly, however, Eq. [5] shows that the 
overall rate of approach to equilibrium will generally be 
dominated by the smallest eigenvalue Xj , independent of 
the initial conditions. Conversely, the largest eigenvalue 
sets the minimum timescale for any significant change in 
the ion population. 



3. RESULTS 

Combining the result from [J2] with existing ioniza- 
tion and recombination rates, we can calculate the ac- 
tual timescales (e-folding times) for an astrophysical 

plasma to reach equilibrium. Figure Q] plots max{l/A}, 
in units of cm -3 s as a function of electron tempera- 
ture for a range of astrophysically abundant elements. 

The companion Figure [2] shows min{l/A}, effectively the 
timescale for any significant change in the ion popula- 
tion. These fi gures use the ionizati on and recombina- 
tion rate from Ma zzotta et al.l (119981). Thes e rates have 
been updated recently (jBrvans et al.ll2009D ; while im- 
portant for spectral studies, the new rates are generally 
within 50% of the older values, not enough to signifi- 
cantly change these results. Figure [T] does show that 
the typical expectation that a plasma will reach equilib- 
rium when n e t ~fewxl0 13 cm -3 s is overly conservative. 
For most atoms, equilibrium will be achieved an order of 
magnitude earlier, or in the case of light ions such as C 
or O at high temperatures, two orders of magnitude or 
more. 

4. DISCUSSION AND CONCLUSIONS 

The complex structure that can be seen in Fig- 
ure [T] is the result of different ions and different ion- 
ization/recombination rates dominating the process at 
different temperatures. An examination of the individ- 
ual ion populations as a function of time showed that the 
most populous ion is often the slowest to reach equilib- 
rium. At lower temperatures (T = f0 4 — I0 6 K), a range 
of L-shell ions, with a range of ionization and recombina- 
tion rates, dominate the final equilibrium state, leading 
to complex variations in the eigenvalues. However, at 
higher temperatures, the curves simplify as the ion pop- 
ulation is dominated by He-like, hydrogenic, or bare ions. 
Figure [5] shows a simpler structure, albeit with one sig- 
nificant 'kink' for each atom. This feature arises as the 
ion population begins to contain a significant number of 
non-neutrals, and is correlated with the first ionization 
potential of each element. 

Although Figure [T] cannot be used to calculate the ion 
population directly, it fills a long-standing need in the X- 
ray spectral community. It is particularly useful in that 
it is entirely independent of initial conditions, equally 
applicable to a neutral medium that is rapidly ionizing 
or tracking the recombination of a fully-stripped plasma. 
Figure [T] can easily be used to determine if the best- 
fit temperature and density-weighted time returned from 
the spectral fitting program implies a significantly non- 
equilibrium plasma or rather one near equilibrium, and 
to see which elements would be the most time-sensitive 
at a given temperature. Figure [U in turn, shows that no 
ion population will be impacted by timescales less than 
~ I0 7 cm~ 3 s and in some cases significantly longer. Al- 
though only collisional processes were included here, the 
method is general. A similar figure could be generated 
for a photoionized plasma held at constant photoioniza- 
tion parameter £, for example when studying a plasma 
near a quiescent X-ray binary after an outburst. 

We thank the anonymous referee for suggesting deter- 
mining the minimum timescale as well as the maximum, 
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Fig. 2. — [Left axis] Density-weighted timcscales (in units of cm 3 s) for C, N, O, Ne, Mg, Al, S, Si, Ar, Ca, Fe, and Ni to achieve one 
c- folding (e~ 1 ) of any change in ionization state for an atom in a constant temperature plasma. [Right axis] Density- weighted timescale 
for an ionization state change of 90%. 
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